California  Institute  of  Technology 


Active  Control  of  Instabilities 
in  Jet  Engines 

Final  Report  for  ONR  N0001492J1677 


John  C.  Doyle 
Caltech  116-81 
Pasadena,  CA  91125 
(818)  395-4808 


19950925  037 


^  DTIC  QUALITY  IIJSPSCTSD  1 

94  8  22  125 


Active  Control  of  Instabilities  in  Jet  Engines 

ONR  N0001492J1677 


Final  Report 
John  Doyle 


1  Summary 

Advances  in  several  areas  of  robust  control  theory  and  its  applications  were  made  during  the 
extension  of  this  program.  Our  emphasis  has  been  in  the  development  of  computable  measures 
of  performance  robustness  for  both  linear  and  nonlinear  systems,  and  in  the  development 
of  computationally  sound  identification  algorithms.  In  the  following  paragraphs  we  briefly 
discuss  the  main  areas  of  our  program.  A  detailed  description  of  the  results  obtained  and 
the  future  directions  to  be  explored  follows.  In  the  detailed  description,  we  will  point  to  the 
different  publications  that  resulted  from  this  research. 

Robust  stability  and  performance  analysis  with  real  parametric  uncertainty  can  be  natu¬ 
rally  formulated  as  a  Structured  Singular  Value,  or  /i,  problem,  where  the  block  structured 
uncertainty  description  is  allowed  to  contain  both  real  and  complex  blocks.  It  is  now  well 
known  that  computation  for  the  general  mixed  ^  problem  is  NP  complete.  Thus,  to  obtain 
acceptable  computation,  we  do  not  attempt  to  solve  the  mixed  /j,  problem  exactly  but  rather 
to  obtain  good  bounds.  The  key  to  obtaining  a  lower  bound  lies  in  the  fact  that  the  n  problem 
may  be  reformulated  as  a  real  eigenvalue  maximization  since  for  any  Q  e  Q,  PriQM)  <  piM). 
The  computational  complexity  of  this  problem  manifests  itself  in  the  fact  that  this  function  is 
non-convex  and  so  it  is  difficult  to  find  the  global  maximum.  Any  local  maximum,  however, 
is  a  lower  bound  to  the  global  maximum.  We  are  currently  working  on  an  efficient  way  to 
compute  a  local  maximum  of  the  this  function  using  a  simple  power  iteration.  We  describe 
this  algorithm  in  more  detail  in  Section  3. 

For  LTI  systems,  performance  tests  are  normally  posed  and  answered  in  the  frequency 
domain,  or  equivalently  on  an  infinite-time  horizon.  Although  this  is  something  of  an  artifice 
for  LTI  systems,  the  nature  and  computation  of  the  tests  are  simplified  as  a  result.  However, 
when  considering  non-linear  or  time- varying  systems  it  is  both  more  natural  and  computa¬ 
tionally  convenient  to  set  the  performance  specification  and  uncertainty  descriptions  in  the 
time  domain,  over  a  finite  horizon.  By  writing  these  requirements  as  quadratic  constraints  we 
can  answer  the  analysis  question:  Is  the  performance  condition  met  for  all  the  plants  in  the 
uncertainty  set?,  by  computing  /r  of  a  constant  matrix,  just  as  in  the  LTI  infinite  horizon  case. 
The  size  of  this  p  test  is  proportional  to  the  length  of  the  horizon,  making  the  computational 
complexity  of  the  lower  and  upper  bound  algorithms  an  important  issue  in  finite  horizon 
tests.  We  have  exploited  the  special  structure  of  the  matrix  associated  with  the  time  domain 
test  in  order  to  modify  the  standard  algorithms  for  the  bounds  and  reduce  the  growth  rate 
of  computation  time  with  the  horizon  length.  This  area  of  our  research  program  is  described 
in  Section  4. 
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Experimental  work  was  also  carried  out  during  the  development  of  this  project.  Advanced 
control  design  techniques  were  applied  to  multimode  stabilization  of  thermally  induced  oscil¬ 
lations  inside  a  cavity,  through  this  experimental  work  we  tested  the  relevance  of  the  different 
theoretical  tools  at  our  reach  to  practical  problems 


2  Preliminaries  and  Notation 

The  notation  we  use  is  as  follows.  7^2,  '^oo  denote  the  Hardy  spaces  of  possibly  vector-  or 
matrix- valued  functions  with  analytic  continuation  on  the  unit  disc,  and  To,  Toe  the  corre¬ 
sponding  Lebesgue  spaces  of  functions  square  integrable  and  essentially  bounded,  respectively, 
on  the  unit  circle,  each  with  norms  |1  •  II2,  ||  •  lU-  T^'^oo  and  TlCoo  are  the  subspaces  of 
and  Too  whose  elements  are  rational  functions.  We  represent  the  integers  by  Z,  the  time  shift 
operator  by  z-S  and  the  identity  matrix  by  I,  where  the  dimensions  will  be  assumed  to  be 
clear  from  the  context,  or  will  otherwise  be  stated.  The  maximum  singular  value  of  .4  is 
denoted  W{A),  and  A*  denotes  the  complex  conjugate  transpose. 

2.1  Linear  Fractional  Transformations  (LFTs) 

We  use  LFTs  to  represent  uncertain  systems  as  shown  in  Figure  1. 


Figure  1:  MD/Uncertain  System 


The  frequency/uncertainty  structure  A  G  A  we  consider  is 

A  =  {diag  , . . . ,  ^slns]  ■  ■  -^2  -^^2} 

BA  =  {a  e  A  :  <  l} 


We  assume  M  =  I  ^  is  a  given  system  realization  matrix  with  A,  B  and  C  parti-  ^ - 

C  D  yCT _ , 

tioned  compatibly  with  the  block  structure  A.  The  input/output  mapping  for  this  system  is  hTJ 

determined  by  ^  | 
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We  refer  to  such  system  models  by  the  pair  (A,M). 

2.2  The  strucutred  singular  value 

For  any  square  complex  matrix  M  we  denote  the  complex  conjugate  transpose  by  M  .  The 
largest  singular  value  and  the  structured  singular  value  are  denoted  by  W{M)  and  n^{M)  re¬ 
spectively.  The  spectral  radius  is  denoted  p{M)  and  Pr{  M)  —  max{\\\  :  \  is  a  real  eigenvalue  of  M] 
with  Pr{M)  =  0  if  M  has  no  real  eigenvalues.  For  any  complex  vector  x,  then  x*  denotes  the 
complex  conjugate  transpose,  and  |x|  the  Euclidean  norm. 

The  definition  of  p  is  dependent  upon  the  underlying  block  structure  of  the  uncertainties, 
which  is  defined  as  follows.  Given  a  matrix  M  €  and  three  non-negative  integers  m,., 

me,  and  me  with  m  :=  m^  -h  +  me  <  n,  the  block  structure  /C(me,  m^,  me)  is  an  m-tuple 
of  positive  integers 

fC.  =r  (/ci ,  .  .  .  ,  fcmr  ?  ?*•••)  f^niy+mc  ^  ^  *  *  *  ^ 

where  we  require  YT=i  ^  dimensions  are  compatible.  This  now  deter¬ 

mines  the  set  of  allowable  perturbations,  namely  define 

Xfc  =  {A  =  block  diag{S[Ik^^ . .  • , ^  ’  *  *  *  ’  * 

<5[  G  G  C,  Af  G  ^4^ 

Note  that  Xfc  e  and  that  this  block  structure  is  sufficiently  general  to  allow  for 

repeated  real  scalars,  repeated  complex  scalars,  and  full  complex  blocks.  The  purely  complex 
case  corresponds  to  —  0.  Note  also  that  the  full  complex  blocks  need  not  be  square  but 
we  restrict  them  as  such  for  notational  convenience. 

Definition  1  ([14])  The  structured  singular  value^  of  a  matrix  M  G  with 

respect  to  a  block  structure  /C(mr,mc,mc)  i^  defined  as 

IJ,)c{M)  =  (  min  {a(BA)  :  det{I  —  BAM)  =  0}\  (5) 

yBAcA'x:  / 

with  =  0  if  no  BA  G  X^c  solves  detfl  —  BAM)  =  0. 

In  order  to  obtain  a  lower  bound  for  g  we  define  the  following  sets  of  block  diagonal 
matrices  (which  are  also  dependent  on  the  underlying  block  structure). 

Q,c  =  {BA  e  Xa:  :  G  [-1  =  1,  Af'Af  =  (6) 

=  {block  diag{e^^^Di^ ....  Dm^ •  T^mr+i  ^  ^  ^  ’  •  ■  •  ’ 

dmchj  ■■  6i  G  [-|  |],0  <  D,  =  D*  G  C‘^-^^',0  <  d,-  G  R}  (7) 
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3  Improved  n  Lower  Bound  Computation 

Robust  stability  and  performance  analysis  with  real  parametric  uncertainty  can  be  naturally 
formulated  as  a  Structured  Singular  Value,  or  fi,  problem,  where  the  block  structured  uncer¬ 
tainty  description  is  allowed  to  contain  both  real  and  complex  blocks.  For  more  engineering 
motivation  of  the  use  of  these  approaches,  see  [1,  16]  and  the  references  therein. 

It  is  now  well  known  that  computation  for  the  general  mixed  n  problem  is  NP  complete 
[52, 13,  9].  It  is  still  a  fundamental  open  question  in  the  theory  of  computational  complexity  to 
determine  the  exact  consequences  of  a  problem  being  NP  complete,  and  we  refer  the  reader  to 
[21]  for  an  in  depth  treatment  of  the  subject.  However,  it  is  generally  accepted  that  a  problem 
being  NP  complete  means  that  it  cannot  be  computed  in  polynomial  time  in  the  worst  case. 
It  is  important  to  note  that  being  NP  complete  is  a  property  of  the  problem  itself,  not  any 
particular  algorithm.  The  fact  that  the  mixed  n  problem  is  NP  complete  strongly  suggests 
that  given  any  algorithm  to  compute  n,  there  will  be  problems  for  which  the  algorithm  cannot 
find  the  answer  in  polynomial  time.  Roughly  speaking,  this  means  that  mixed  p  cannot  be 
computed  exactly  in  the  worst  case  without  entirely  unacceptable  growth  in  computation  cost 
with  problem  size.  For  all  practical  purposes  even  moderately  large  examples  of  such  problems 
are  computationally  intractable.  To  obtain  acceptable  computation,  we  do  not  attempt  to 
solve  the  mixed  /j,  problem  exactly  but  rather  to  obtain  good  bounds.  Furthermore,  [13] 
suggests  that  even  approximate  methods  are  also  NP  complete,  so  we  will  not  expect  good 
worst  case  behavior  but  rather  aim  for  good  typical  behavior. 

The  key  to  obtaining  a  lower  bound  lies  in  the  fact  that  the  /.i  problem  may  be  reformulated 
as  a  real  eigenvalue  maximization  [66]. 

maxpriQM)  =  p{M)  (8) 

QeQ 

Here  Pr  is  the  real  spectral  radius,  and  Q  is  a  set  of  matrices  with  appropriate  structure. 
This  immediately  gives  us  a  theoretical  lower  bound  since  we  have  that  for  any  Q  G  Q, 
Pr{QM)  <  p{M).  The  computational  complexity  of  this  problem  manifests  itself  in  the  fact 
that  this  function  is  non-convex  and  so  it  is  difficult  to  find  the  global  maximum.  Any  local 
maximum,  however,  is  a  lower  bound  to  the  global  maximum. 

The  idea  then  is  to  find  an  efficient  way  to  compute  a  local  maximum  of  the  function 
PriQM)  over  Q  e  Q.  It  turns  out  that  this  can  be  done  by  means  of  a  simple  power  iteration. 
The  iteration  scheme  usually  converges  fairly  rapidly  and  each  iteration  of  the  scheme  is 
very  cheap,  requiring  only  such  operations  as  matrix-vector  multiplications  and  vector  inner 
products. 

Unfortunately  the  lower  bound  power  iteration  is  not  always  guaranteed  to  converge. 
Although  one  can  still  obtain  a  lower  bound  in  this  case,  it  may  no  longer  correspond  to  a 
local  maximum  of  the  real  spectral  radius,  and  so  the  bound  may  be  poor.  While  Branch  and 
Bound  techniques  are  an  effective  way  of  obtaining  improved  bounds,  results  in  [41]  strongly 
suggest  that  the  quality  of  the  upper  and  lower  bounds  for  p  are  critical  to  the  performance 
of  any  Branch  and  Bound  scheme  to  refine  them.  This  means  that  any  improvement  in  the 
lower  (or  upper)  bound  performance  is  highly  desirable,  whether  one  wishes  to  use  the  bounds 
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directly,  or  as  part  of  a  Branch  and  Bound  scheme. 

There  is  additional  motivation  for  further  development  of  the  lower  bound  computation: 
there  are  important  problems  for  which  computation  will  require  generalizations  of  the  current 
lower  bound  algorithm  and  these  generalizations  are  expected  to  be  more  difficult  in  the 
sense  that  it  will  be  harder  to  converge  to  local  maximums.  One  such  problem  is  the  Model 
Validation  problem,  where  we  wish  to  determine  whether  measured  data  (time  histories) 
are  consistent  with  a  given  robust  control  model  with  pre  specified  bounds  on  uncertainty 
operators  and  unknown  noise.  Another  is  the  robust  control  ID  problem.  Here,  we  wish  to 
determine  parameters  that  best  fit  measured  data  in  that  they  select  a  robust  control  model 
with  (in  some  sense)  minimal  uncertainty  or  noise.  Both  of  these  problems  are  distinguished 
by  the  presence  of  known  signals.  Problems  of  this  nature  in  general  require  an  extension  of 
the  current  lower  bound  computation.  Note  that  the  lower  bound  algorithm  finds  solutions 
to  the  associated  loop  equations  and  for  the  Model  Validation  problem  and  the  ID  problem 
these  solutions  include  noise  signals  and/or  parameters  of  the  model. 

The  increased  difficulties  encountered  in  the  mixed  ^  problem  stem  from  a  simple  fact.  In 
the  complex  only  case,  the  real  spectral  radius  maximization  is  equivalent  to  a  spectral  radius 
maximization.  That  is,  if  we  maximize  the  spectral  radius,  then  we  can  easily  create  a  real 
eigenvalue  equal  to  the  spectral  radius.  Hence  we  may  use  an  algorithm  that  generalizes  the 
Rayleigh-Ritz  method  of  finding  the  spectral  radius.  This  is  not  the  case  in  the  general  mixed 
p.  problem.  Here  we  wish  to  maximize  the  largest  real  eigenvalue  even  when  the  spectral 
radius  is  much  larger.  Consequently  local  maxima  are  not  generally  stable  equilibria  of  the 
standard  power  iteration.  This  is  because  errors  with  a  component  in  the  direction  of  the 
eigenvector  corresponding  to  the  eigenvalue  that  achieves  the  spectral  radius  are  amplified.  A 
consequence  of  this  is  that  the  lower  bound  algorithm,  as  implemented  in  the  /i-tools  toolbox, 
fails  to  find  a  local  maximum  on  a  significant  number  of  problems. 

In  [61]  it  is  demonstrated  that  there  is  potential  for  substantial  improvement  over  the 
standard  algorithm.  The  present  goal,  and  the  focus  of  the  lower  bound  work  over  the  past 
year,  is  to  develop  an  efficient  algorithm  that  retains  the  performance  of  the  current  algorithm 
when  this  is  satisfactory,  and,  when  is  it  not,  find  a  (large)  local  maximum  as  efficiently  as 
possible. 

One  aspect  of  this  research  is  the  investigation  and  development  of  a  shift  and  inverse 
algorithm.  This  is  a  generalization  of  shift  and  inverse  algorithms  for  finding  eigenvalues  and 
eigenvectors.  Observe  that  if  Xx  —  Mx,  then  (A  —  f3)x  =  [Ad  —  (3I)x  and  (A  —  /3)  x  — 
[Ad  -  (3d)-^x.  If  /?  is  close  enough  to  A,  then  [(A  -  /?)"^|  is  equal  to  the  spectral  radius  of 
[Ad  -  Pdy^  and  a  power  algorithm  would  tend  to  find  the  eigenvalue  (A  -  ^)“^  which  is 
equivalent  to  finding  the  eigenvalue  X  of  Ad . 

This  technique  has  been  used  for  the  mixed  jj,  lower  bound  real  eigenvalue  maximization 
problem.  The  idea  is  to  chase  an  eigenvalue  to  a  local  maximum:  we  shift- and-inverse  iterate 
to  find  an  eigenvalue  of  MA,  then  we  update  A  (this  is  moving  in  the  domain  of  optimization  ) 
based  on  the  eigenvectors  and  repeat.  The  difficulty  with  this  approach  is  that  A  must  change 
slowly  so  that  the  next  shift  is  close  to  the  eigenvalue  we  are  chasing.  This  necessitates  many 
steps  and  the  steps  are  not  as  cheap  as  we  would  like  (the  cost  grows  with  the  cube  of  the 
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size  of  M),  so  this  method  is  less  than  ideal. 

Another  approach  is  based  on  the  fact  that  the  optimum  over  the  parameters  jointly  is 
also  an  optimum  over  a  subset  of  the  parameters  when  the  complement  of  the  subset  is  held 
fixed  at  their  optimal  values.  Now  the  optimum  over  the  complex  variables  is  quick  and  easy: 
this  is  the  standard  power  algorithm  on  a  strictly  complex  problem.  The  optimization  over 
the  real  parameters  may  be  approached  with  conventional  optimization  techniques  as  all  the 
real  blocks  are  scalar  times  identity  blocks  and  thus  there  are  relatively  few  real  parameters. 
The  drawback  to  this  problem  is  that  the  conversion  between  the  two  optimization  problems 
is  not  cheap,  and  conventional  optimization  techniques  still  may  require  many  steps. 

While  the  current  implementation  of  these  ideas  results  in  an  algorithm  that  is  much  better 
than  the  standard  algorithm  in  that  it  finds  much  better  solutions  to  most  of  the  difficult 
problems  reasonably  quickly,  we  believe  there  is  room  for  further  substantial  improvement. 
These  improvements  would  come  not  only  from  refinements  to  the  algorithms  described  above, 
but  also  from  a  careful  combination  of  these  algorithms  and  the  standard  power  algorithm. 


4  Finite  Time  Horizon  Robust  Performance  Analysis 

4.1  Introduction 

For  LTI  systems,  it  is  convenient  to  pose  performance  tests  over  an  infinite  horizon.  This 
artifice,  in  this  particular  subset  of  systems,  simplifies  rather  than  complicates  the  nature  and 
computation  of  the  tests.  However,  if  we  want  to  work  with  non-linear  or  time- varying  systems 
it  is  more  natural,  and  computationally  convenient,  to  set  the  performance  specification  and 
uncertainty  descriptions  in  the  time  domain,  over  a  finite  horizon.  As  in  the  LTI  infinite 
horizon  case,  when  these  requirements  are  written  as  quadratic  constraints  we  can  answer  the 
analysis  question:  Is  the  performance  condition  met  for  all  the  plants  in  the  uncertainty  set?, 
by  computing  the  structured  singular  value  /x  of  a  constant  matrix.  The  size  of  the  /x  test 
required  grows  linearly  with  the  length  of  the  horizon,  making  the  computational  complexity 
of  the  lower  and  upper  bound  algorithms  an  important  issue  in  finite  horizon  tests.  By 
exploiting  the  special  structure  of  the  matrix  associated  with  the  time  domain  test,  we  can 
modify  the  standard  algorithms  for  the  bounds,  to  reduce  the  growth  rate  of  computation 
time  with  the  horizon  length. 

When  the  finite  time  tests  are  applied  to  LTI  systems  over  increasing  horizons,  we  would 
expect  to  recover  the  infinite  horizon  results.  The  connection  between  the  exact  tests  is  a 
previously  established  property  of  the  structured  singular  value.  By  using  losslessness  of  the  S- 
Procedure  theorem  (see  [37]  and  the  references  therein)  we  have  extended  these  connections  to 
the  upper  bounds  of  the  finite  time  and  frequency  domain  tests.  We  expect  these  connections 

to  shed  new  light  on  the  nature  of  both  tests. 

One  of  the  most  exciting  applications  of  this  theory  is  the  analysis  of  non-linear  systems. 
We  can  use  the  analysis  tool  for  finite  horizon  linear  time  varying  (FHLTV)  systems  as  one 
step  in  an  iterative  procedure  for  analysis  of  non-linear  systems  along  trajectories.  This 
procedure,  which  alternates  simulation  of  the  nonlinear  system,  and  linearization  around  a 
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perturbed  trajectory  can  locally  find  worst  case  noise  signals.  Ongoing  research  is  dedicated  to 
the  development  of  efficient  and  numerically  sound  algorithms  for  carrying  out  this  procedure. 

4.2  Problem  Setup 

4.2.1  Linear  Time  Varying  Systems  in  Finite  Horizon 

LTI  systems  are  normally  described  as  transfer  functions,  namely  a  map  from  U  into  U. 
However,  in  order  to  derive  computable  tests  we  have  to  describe  them  in  terms  of  constant 
matrices.  This  can  be  done  either  in  state  space,  by  incorporating  the  delay  operator  in  the 
uncertainty,  or  in  the  frequency  domain,  by  doing  a  search  over  frequency.  The  analysis  at 
each  frequency  point  reduces  to  a  constant  matrix  problem.  None  of  these  approaches  extends 
naturally  to  time  varying  systems  considered  over  a  finite  time  horizon.  However,  a  FHLTV 
svstem  can  be  represented  as  a  finite  dimensional  map  (i.e.  a  constant  matrix)  by  mapping 
the  time  axis  into  a  spatial  axis.  To  illustrate  this  point  let’s  consider  a  discrete  time  linear 
time  varying  system,  over  a  time  horizon  of  k  steps.  The  system  is  described  by  the  equations. 

Xi+i  =  A(i)xi  +  B{i)ui  (9) 

y^  =  C(i)xi  +  D(i)ui  i  =  0-  ■  -k  -  1 

These  equations  can  be  rewritten  as: 

-)  =  M(.)(:;)  .  =  (10) 

We  can  define  a  map  from  the  initial  state  and  the  time  history  of  the  inputs  (from  t  =  0 
to  t  =  A:  -  1),  to  the  final  state  and  the  time  history  of  the  outputs  (from  t  -  0  to  t  =  k  -  1). 
We  denote  this  mapping  with  the  symbol  M[k].  (See  Figure  2.)  The  system  is  now  represented 
by  a  single  constant  matrix  M[k],  over  which  we  can  write  our  performance  bounds.  Note  that 
the  dimensions  of  the  matrix  grow  linearly  with  the  length  of  the  time  horizon.  This 
requires  that  a  strong  emphasis  be  put  on  the  development  of  efficient  algorithms  for  the 
computation  of  the  stability  tests. 

4.2.2  Time  Domain  Uncertainty  -  Performance  Specifications 

In  finite  time  horizon  problems,  performance  is  the  only  relevant  concept,  since  stability  is  just 
an  artifact  of  infinite  time.  For  this  work  we  consider  performance  and  uncertainty  conditions 
which  can  be  written  as  quadratic  constraints  in  the  variable: 

2  =  (xo  uo  •  ■  •  Ui-1  .rj.  Vo  ■■■  2/t-i)* 

The  uncertainty  and  performance  constraints  have  the  form: 

(Ti{z)  >  0  (11) 

Where  t  =  0  for  the  performance  constraint,  and  i  =  1 . .  .n  for  the  uncertainty  constraints. 
In  defining  these  QCs  we  use  the  following  conventions: 


( 


Figure  2:  Finite  time  horizon  system  as  a  constant  matrix  transformation 


•  When  the  QC  describes  an  uncertainty  condition,  a, (2)  >  0  implies  that  there  exists  an 
instance  of  the  uncertainty  with  the  desired  properties  that  can  account  for  the  data  in 

z. 

•  When  the  QC  describes  the  performance  specification,  <70(2)  >  0  implies  that  the  data 
does  not  meet  the  performance  required. 

As  seen  in  [60]  and  [48],  there  is  a  A-structure  associated  with  this  uncertainty-performance 
description.  We  give  two  examples  of  how  to  convert  from  one  representation  to  the  other. 
The  performance  quadratic  constraint: 

i:ii<.(ft)p+iwr>i:ii!'Wip+iKii’  (12) 

2=0  *=1 

is  met  if  the  total  energy  at  the  output  is  smaller  than  the  total  energy  at  the  input.  The 
corresponding  A  structure  for  this  QC  is  a  map  from  the  final  state  and  the  full  time  history 
of  the  output  to  the  initial  state  and  the  full  time  history  of  the  input  as  shown  in  Figure  3 
If  the  uncertainty  we  wish  to  describe  is  a  bounded  varying  gain,  we  use  the  QCs: 

|l2/(0ir- IKOir  >  0  0<i</c-l  (13) 

The  associated  uncertainty  structure  is  block  diagonal,  as  shown  in  Figure  4. 

4.2.3  Time  Domain  [a. 

Using  the  setup  described  in  the  preceding  section,  it  is  easy  to  see  that  the  time  domain 
performance  problem  can  be  recast  as  a  standard  jA  problem  on  the  matrix  with  respect 
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Figure  3:  Total  output  energy  to  total  input  energy  performance 
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Figure  4:  Instant  by  instant  output  energy  to  input  energy  uncertainty 


to  a  block  structure  [A^^,  Ap]  that  can  be  derived  from  the  preceding  definitions.  Therefore 
the  performance  specifications  are  met  robustly  whenever  the  following  condition  holds: 

IJ-[A.,A,]{M[k])  <  1 

where  ^  is  the  standard  structured  singular  value  as  defined  for  example  in  [42]. 

Thus  we  have  the  following  definition: 

Definition  2  Time  domain  ^  Given  a  system  of  the  form  (10),  observed  over  a  time  horizon 
of  k  steps,  and  uncertainty  specifications  given  by  the  quadratic  constraints  Oi,  i  ==  0  —  *n 
then: 


4.3  Necessary  and  Sufficient  Conditions 
4.3,1  A  Necessary  condition  for  performance 

If  we  scale  the  variable  z  to  Zx*  i.e., 


z  =  (a:o  uo  •  *  ■  Uk-i  Xyo  •  •  • 

and  write  the  performance-uncertainty  QCs  in  we  will  be  relaxing  the  requirements:  the 
bigger  A  is  the  easier  it  is  for  the  system  to  meet  these  requirements.  Let  A^ax  be  the  maximum 
value  of  A  for  which  the  conditions  are  not  met.  A  necessary  condition  for  performance  is 
Amax  <  1-  As  shown  in  [44]  A^ax  =  va3ixp{QM[k]),  where  Q  is  unitary  and  has  the  same 
structure  as  [A^,  Apj.  We  can  not  in  general  compute  this  value,  but  we  can  compute  local 
maxima  A  of  p{QM[k])  along  with  the  signals  [xq  Uq  •  •  •  for  which  the  conditions 

9 


I 


fail.  This  not  only  gives  us  a  way  to  compute  a  necessary  condition  for  performance  (v.z. 

X  <  1),  but  it  also  gives  us  a  way  of  finding  worst  case  signals  for  the  system.  This  fact  wdll  be 
important  when  applying  this  theory  to  the  analysis  of  nonlinear  systems  (See  Section  4.5). 

Computationally  sound  algorithms  have  been  developed  (see  [66]  and  references  therein) 
to  compute  a  lower  bound  for  the  maximum  over  Q  of  the  spectral  radius  of  QM.  In  order 
to  apply  this  algorithm  to  our  problem  some  modifications  were  made  in  order  to  exploit  the 
structure  of  the  matrix  M[k]  and  avoid  quadratic  growth  in  computation  time  with  the  length 
of  the  horizon  ([60]). 

4.3.2  Sufficient  conditions  for  performance 

As  with  the  necessary  condition,  we  have  derived  a  sufficient  condition  for  robust  performance 
from  the  constant  matrix  fi  problem  described  in  section  4.2.3. 

Definition  3  Time  domain  upper  bound 

tdub{„.]{M{-),k)  =  w6[A„,Ap](M[fc]) 

where  ub  denotes  the  standard  p  upper  bound. 

It  is  immediate  that  tdw6{^.}(M(-),  fe)  <  1  is  a  sufficient  condition  for  robust  performance. 
The  same  computational  complexity  problems  arising  in  the  lower  bound  arise  in  the 
computation  of  the  upper  bound.  However,  it  is  not  as  straightforward  to  use  the  structure 
in  the  matrix  M[k]  to  reduce  the  growth  of  the  computation  time,  when  using  current  state 
of  the  art  optimization  algorithms  for  LMIs. 

Using  gradient  search,  we  can  develop  an  algorithm  that  is  slower  on  small  problems  when 
compared  to  the  LMI  methods,  but  whose  computation  time  doesn’t  degrade  as  much  with 
the  number  of  time  samples.  However  our  tests  of  this  algorithm  are  still  preliminary  and 
more  extensive  experimentation  and  further  research  is  being  completed  in  this  area. 

4.4  Connections  to  the  frequency  domain  tests 

Assuming  the  system  under  consideration  is  LTI,  the  uncertainty  description  is  repeated  at 
each  time  instant  and  the  performance  condition  is  given  as  a  full  block  mapping  the  final 
state  to  the  initial  state,  we  have  developed  connections  between  the  time  domain  tests  and 
the  corresponding  frequency  domain  tests.  The  following  two  theorems  give  an  idea  of  the 
nature  of  this  relationship. 

Theorem  1  [44]  Given  a  system  M  an  uncertainty  structure  A^,  and  e  >  0,  then  there  exists 
A'  >  0  such  that  for  all  k  >  K,  tdp{M{-),  k)<l  +  e,  if  and  only  if  Pa{M)  <  1 

When  the  system  is  a  constant  matrix,  we  use  results  on  the  losslessness  of  the  S-procedure 
([37]  ,  [6],  [48])  to  generalize  the  preceding  result  to  the  upper  bounds  in  the  frequency  and 
time  domains,  respectively: 
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Theorem  2  Given  a  system  M  an  uncertainty  structure  set  A^,,  and  e  >  0,  then  there  exists 
A'  >  0  such  that  for  all  k  >  K,  tdub(M{-),  k)  <  I  +  e  if  only  if  fdab^^iM)  <  1 

From  these  connections  we  expect  to  derive  a  better  understanding  of  the  nature  of  the 
frequency  domain  tests  for  LTI  systems,  especially  in  view  of  the  results  in  [49]  and  [50]. 


4.5  Application  to  Non-Linear  Systems 

An  exciting  application  of  this  theory  is  as  part  of  an  iterative  algorithm  to  analyze  the  per¬ 
formance  of  discrete  time  non-linear  systems.  The  state  of  the  art  for  non-linear  performance 
analysis  is  extensive  simulation  under  different  perturbations  and  uncertainty  instances.  We 
propose  an  iterative  algorithm  that  can  locally  find  worst  case  perturbations  and  uncertainty. 
This  algorithm  can  be  used  to  reduce  the  number  of  simulations  needed  in  the  analysis  of  the 
plant. 

4.5.1  An  iterative  analysis  tool 

We  consider  a  discrete  time  non-linear  system  described  by  the  equations: 

Xi+i  =  A{xi)xi  +  Bi{xi)ui  +  B2ixi)wi  (14) 

yi  =  C{xi)xi  +  Di{xi)ui  +  D2{xi)wi  i  =  l...k 

We  analyze  this  system  along  the  trajectory  generated  by  a  previously  specified  input  u. 
Our  objective  is  to  find  the  norm  bounded  disturbance  signal  w  that  locally  maximizes  the 
distance  from  the  perturbed  to  the  nominal  trajectory  measured  in  the  2-norm. 


Nominal  Trayectory 


Perturbed  Trayectory 
After  one  iteration. 


Figure  5:  First  step  in  Simulation 


To  start  the  procedure  we  simulate  the  nonlinear  system  along  the  nominal  input,  in  the 
absence  of  noise.  We  then  consider  perturbations  of  the  system  along  this  trajectory  caused 
by  the  disturbance  w. 


+  «Si  =  +  Bi{xi)ui  -f  Boixifwi  (15) 

j/W  +  =  C(xi)xi  -f  Diixi)ui  +  D2ixi)wi  i=l---k 
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where  we  have  decomposed  a:  and  y  into  a  nominal  and  perturbed  component.  The  pertur¬ 
bations  around  the  nominal  will  verify  the  following  uncertain  time  varying  linear  equations: 


/  \ 


\  / 


Wi  \ 
Wk  / 


(16) 


These  equations  have  the  form  described  in  the  proceeding  sections.  We  can  thus  use  the 
td^i  lower  bound  algorithms  to  obtain  worst  case  perturbations  for  this  system.  We  scale 
this  perturbation  to  the  predetermined  noise  size,  and  simulate  the  original  non-linear  sys¬ 
tem.  Denote  by  wq  the  time  history  of  the  noise  signal,  w,  and  i/q  the  time  history  of  the 
perturbation  around  the  nominal  trajectory. 

We  can  repeat  the  previous  linearization  around  the  perturbed  trajectory.  If  we  denote 
by  w  -t-  luo  and  y^^  -t-  y^^  the  total  noise  and  total  distance  to  the  nominal,  the  equations  for 
the  new  linearization  are:  ,  _  . 


A  * 


-Pll  ^12 
P2I  P22 


W 


(17) 


.  •  1  •  II  1*1 

Our  objective  is  now  to  find  the  perturbation  w  that  maximizes  the  gain  ||to+iDo||  ’  ^ 

fixing  the  size  of  the  total  noise  w  +  Wo  and  of  the  uncertainty.  Although  this  appears  to  be 
of  a  somewhat  different  nature  than  the  test  in  16,  we  can  rewrite  it  to  put  it  in  the  same 
setup  as  before.  This  is  made  clear  by  the  block  diagram  in  Figure  6. 


Figure  6:  Setup  for  the  second  step  in  the  simulation 


In  this  new  problem  the  sizes  of  A„  and  A„  are  fixed  and  we  are  trying  to  minimize  the 
size  of  Ap. 

After  the  second  step  we  simulate  the  nonlinear  system  with  the  new  perturbation.  I'Ve 
check  at  this  point  if  the  linearization  in  the  previous  step  is  still  valid  in  the  new  trajectory. 
If  it  is,  we  stop  the  algorithm;  otherwise  we  repeat  the  second  step. 
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4.5.2  Computational  Aspects 

Although  most  of  the  theoretical  issues  in  the  development  of  this  algorithm  have  been  com¬ 
pleted,  more  research  is  needed  to  fully  understand  the  behavior  of  the  practical  algorithms. 
The  algorithms  are  of  similar  nature  to  the  ones  used  in  the  usual  performance  analysis  tests. 
However  they  are  applied  to  a  different  set  of  matrices.  We  still  do  not  have  enough  knowledge 
of  the  convergence  characterisitcs  of  the  algorithms  over  this  new  set.  Our  current  efforts  are 
in  this  direction. 

5  Model  Validation  and  System  Identification 

5.1  Introduction 

The  basic  element  of  any  quantitative  approach  to  scientific  and  technological  problems  is  a 
mathematical  model.  These  models  are  typically  obtained  using  some  combination  of  first 
principles  and  identification  from  experimental  data.  First  principles  models  build  on  scientific 
knowledge  which  is  itself  often  the  consequence  of  extensive  experimentation.  The  most  useful 
models  would  in  principle  arise  from  a  judicious  combination  of  first  principles  and  system 
identification,  with  the  measure  of  a  models  usefulness  being  highly  context  dependent.  In 
general,  we  might  start  with  a  model  with  some  a  priori  structure  (here  we  could  use  some 
first  principles  knowledge),  and  which  includes  some  description  of  the  a  priori  uncertainty 
(parameters,  disturbances,  etc.).  After  performing  an  experiment,  we  pose  the  mathematical 
problem  of  finding  values  of  the  uncertainty  that  fit  our  data. 

An  extensive  field  of  research  has  been  built  in  pursuing  the  answer  to  this  problem;  a 
standard  reference  is  [28].  A  canonical  example  of  the  methods  of  standard  system  identifi¬ 
cation,  is  the  fitting  of  a  parametric  model  by  using  prediction  error  methods  (PEMs).  We 
assume  the  following  model  structure, 

y  ^  G(\,e)u+  H{X,9)n  (18) 

where  A  is  the  shift  operator,  ^  is  a  vector  of  parameters,  G  and  H  are  discrete  time  systems, 
and  n  is  assumed  to  be  noise.  Given  data  for  u,y,  these  methods  attempt  to  find  values  of 
9.  n  which  explain  the  data;  since  many  solutions  typically  exist,  the  standard  approach  is  to 
choose  the  solution  which  minimizes  some  norm  of  n. 

A  related  problem  is  the  model  validation  question:  given  a  model  and  data,  does  the 
model  explain  the  data?  In  the  previous  example,  the  typical  situation  is  that  values  of  9 
have  already  been  chosen,  and  we  wish  to  determine  whether  the  model  explains  a  set  of 
data  with  a  “plausible”  instance  of  n.  Although  these  two  questions  appear  to  be  different 
from  the  point  of  view  of  classical  system  identification,  if  we  consider  the  various  forms  of 
uncertainty  in  an  equal  footing,  this  distinction  disappears,  as  will  now  be  explained  .  The 
basic  mathematical  question  that  needs  to  be  answered  is  the  following: 

Q:  Given  an  uncertain  model  and  experimental  data,  do  there  exist  values  of  the 
uncertainty  such  that  the  model  fits  the  data? 
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In  the  above  example  of  system  identification,  the  uncertainty  used  to  explain  the  data 
is  made  up  of  real  parameters  and  a  disturbance,  whereas  in  the  model  validation  example 
only  the  disturbance  is  used;  but  the  essence  is  the  same,  and  one  could  think  of  including 
other  descriptions  of  uncertainty  (e.g.  unmodeled  dynamics)  into  the  problem.  In  PEM,  the 
only  feature  which  is  added  to  the  previous  question,  is  an  optimization;  among  the  possible 
combinations  of  parameters  and  noise  that  give  an  affirmative  answer  to  Q,  the  minimization 
of  the  noise  is  desired.  Also,  parameters  and  noise  are  treated  in  the  PEMs  in  a  non-symmetric 
way:  the  parameter  values  obtained  will  be  considered  fixed  in  the  resulting  model,  but  not 
the  actual  values  of  the  noise. 

While  this  may  be  a  reasonable  procedure  under  many  circumstances,  it  is  not  objectively 
clear  that  it  is  the  best  approach:  conceivably,  one  could  prefer  a  final  model  where  some 
parametric  uncertainty  is  left.  These  choices  are  up  to  the  user,  and  more  choices  will  arise 
as  we  include  other  sources  of  uncertainty.  Also,  we  can  imagine  multiple  experiments,  where 
some  of  the  uncertainty  is  fixed  to  have  a  common  value  across  experiments,  whereas  other 
parts  (e.g.,  noise,  parametric  variations  due  to  changes  in  experimental  conditions,  unmodeled 
dynamics),  is  allowed  to  vary  from  one  experiment  to  another. 

Therefore,  a  general  methodology  for  model  validation  and  system  identification  should 
provide  computational  tools  for  answering  the  above  general  question  Q  for  very  rich  uncer¬ 
tainty  structures,  including  noise,  unmodeled  dynamics,  and  parameters.  Even  the  parameters 
may  have  a  variety  of  allowable  variations  across  experiments  or  across  time  within  one  ex¬ 
periment.  On  top  of  this  basic  tool,  one  can  think  of  superimposing  other  criteria  to  select  the 
final  model,  or  to  perform  the  validation  question  iteratively  on  various  structures  and  sizes 
of  parameters  to  select  the  final  model.  We  will  concentrate  in  what  follows  on  mathematical 
and  computational  machinery  to  address  Q. 

5.2  Model  validation/ID  in  an  LFT  setting 

5.2.1  An  input-output  approach 

IVhen  the  uncertainty  entering  a  model  is  described  in  an  LFT  manner,  techniques  for  ro 
bustness  analysis  can  be  adapted  to  answer  the  general  model  validation  question  Q.  We 
begin  by  reviewing  briefly  the  approach  given  in  [40],  to  address  this  question,  which  involves 
a  generalization  of  the  structured  singular  value  n- 

A  model  with  incorporated  uncertainty  is  illustrated  in  Fig.  7.  For  simplicity,  there  are 
no  dynamics  in  this  model,  all  the  signals  in  the  diagram  are  constant  vectors,  and  F,  A 
constant  matrices.  We  can  think  of  identifying  a  static  map;  also,  in  the  case  of  system 
identification,  there  are  ways  of  rewriting  the  equations  with  no  explicit  time  horizon;  this 
will  be  mentioned  later. 

The  model,  P,  is  given,  as  is  the  perturbation  structure  for  A.  The  signals  d  and  n  rep¬ 
resent  unknown  disturbances  and  noise  respectively.  The  model  carries  with  it,  assumptions 
on  the  size  of  the  unknown  components.  In  particular  it  is  assumed  that  A  G  ||d||  ^  1, 

and  ||u||  <  1. 
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Figure  7:  Block  diagram  of  the  model  validation  problem 


The  LFT  form  of  the  input-output  description  of  this  system  is, 

y  +  Wnn  =  A*p[^ 

The  perturbation,  A,  has  a  prescribed  block  structure  and  norm  bound:  A  G  BA.  We  will 
assume  that  the  models  under  consideration  are  well  posed  for  all  A  6  BA.  Equivalently, 
yu( Pii)  <  1.  This  ensures  that  the  inverse  in  the  fractional  equations  is  always  well  defined. 

An  input-output  experiment  has  been  run  and  a  measured  datum  {y,u)  is  available.  The 
model  validation  problem  is  to  determine  whether  there  is  an  element  of  the  model  set,  which 
includes  unknown  perturbations  and  unknown  noise  and  disturbance  signals,  such  that  the 
observed  datum  can  be  produced  by  the  model.  This  is  a  necessary  condition  for  the  model 
to  be  able  to  describe  the  system.  In  the  work  presented  here  we  will  consider  only  the 
case  where  P  and  A  are  constant,  complex  valued  matrices.  Application  to  the  more  useful 
problem  involving  dynamic  systems  is  given  in  [57].  The  problem  can  be  stated  as  follows. 

Problem  3  (Model  Validation)  Assume  that  m(Ai)  <  1-  Does  there  exist  a  A  6  BA, 
and  signals  d  and  n,  satisfying  ||d|l  <  1  and  Hull  <  1,  such  that 

y  +  W„n^A^p[i] 

Any  {A,d,n)  satisfying  these  conditions  will  be  referred  to  as  admissible.  We  will  assume 
that  Wn  is  invertible.  This  means  that  every  measured  output,  y,  is  modeled  as  having  some 
non-zero  noise  added  to  it.  Again  this  is  reasonable  in  any  physically  motivated  problem. 
Note  that  if  y  -  P23U  =  0  then  the  model  validation  problem  is  solved  trivially.  Assuming 
that  this  is  not  the  case,  the  problem  can  be  reduced  to  a  generalized  y  test: 
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Theorem  4  There  exists  a  A  e  BA,  and  signals  d  and  n,  with  PH  <  1  and  ||n||  <  1 
satisfying 

y  +  W„n  =  F,{P,A)[i 

iff^l,iP)>l. 

In  the  theorem, 

Pll  Pl2  Pl3'*i 

P=  0  0  1 

_W-^P21  W-^P22  W-\P2sU-y)_ 

is  a  matrix  obtained  from  the  model  and  the  data,  and  /^^(•)  is  a  matrix  function  which 
generalizes  //.  ^Ve  refer  to  [40]  for  the  definition  of  this  function  and  details  of  the  proof  of 
these  results. 


5.2.2  Model  Validation  as  Behavioral  analysis 

A  remarkable  recent  development  [15]  is  that  the  behavioral  robustness  analysis  framework 
provides  an  alternative  (equivalent)  way  of  posing  this  model  validation  question;  since  this 
framework  is  more  closely  related  to  the  rest  of  our  program  we  will  concentrate  in  what 
follows  on  this  new  formulation. 

Figure  8  shows  a  generic  input-output  Model  Validation  structure,  slightly  more  general 
than  Figure  7.  z  is  a  vector  of  unknown  inputs,  constrained  by  ||2||  <  1,  which  could  encom¬ 
pass  d  and  n  in  the  previous  structure.  A  varies  in  BA,  and  we  assume  //(Pn)  <  1  as  before, 
u,  y  are  the  measured  inputs  and  outputs,  (we  assume  that  y  ^  P22w)-  The  model  validation 
question  is,  once  more,  to  find  values  of  (A,  2)  which  account  for  u,y. 


Figure  8:  A  standard  input-output  model  validation  setup. 


In  figure  9,  the  same  equations  are  represented  in  output  nulling  form  (see  [12]).  Com¬ 
bining  u  and  y  into  a  vector  w  (which  includes  all  the  observed  data)  gives  the  behavioral 
representation  of  figure  10.  Note  that  the  input-output  partition  has  been  eliminated  from 
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the  model;  in  fact,  it  could  very  well  be  that  we  wish  to  validate  some  model  based  on  ob¬ 
servations  of  a  system  where  this  distinction  is  not  available  to  us,  and  we  use  a  behavioral 
model  as  a  starting  point.  We  can  now  incorporate  the  data  into  the  matrix  in  a  similar 
manner  to  what  was  done  in  the  input  output  setting,  leading  to  figure  11.  Up  to  this  point, 
no  special  advantage  was  obtained  by  using  the  behavioral  equations. 


Figure  9:  Reduction  to  behavioral  form 


0 


w 

z 


Figure  10:  A  standard  behavioral  model  validation  setup. 


We  have  an  analysis  problem  involving  a  constant  matrix,  but  we  must  deal  with  the 
norm  constraint  on  the  disturbance  2.  At  this  point  the  behavioral  framework  gives  us  the 
possibility  of  expressing  this  constraint  in  GON  form,  in  the  same  way  as  in  the  case  of 
IQCs  (see  [48]).  We  enforce  the  constraint  \\z\\  <  1  by  adding  two  equations  to  the  diagram. 
(See  Figure  12.)  The  first,  p  =  Apl,  introduces  a  variable,  p,  used  on  the  output  side  of  a 
norm  bounded  block.  The  second  ensures  that  2  =  p.  With  ||Ap||  <  1,  we  now  have  that 
||2:||  <  1.  In  Figure  12  the  global  GON  representation  for  the  model  validation  problem  is 
depicted.  The  fact  that  one  signal  is  constrained  to  be  1  is  irrelevant,  since  everything  can 
be  normalized  by  linearity.  Therefore  the  model  validation  question  reduces  to  the  analysis 
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Figure  11:  The  model  validation  setup  with  no  signal  constraints. 


question  ‘‘do  there  exist  nontrivial  signals  satisfying  these  GON  equations  for  A  G  BA  ?  , 
which  can  be  addressed  using  the  machinery  of  [48].  In  this  case  the  D  matrix  of  the  GON 
representation  is  tall,  and  will  be  full  column  rank  unless  P22W  =  0,  which  is  a  trivial  case 
because  the  nominal  model  satisfies  the  data.  Therefore  it  can  be  reduced  to  the  question 


Figure  12:  Model  validation  as  an  analysis  problem 


The  LMI  developed  in  [48]  will  be  a  sufficient  condition  for  invalidating  the  model:  no 
signals  exist  and  therefore  the  data  cannot  be  explained. 

Necessary  conditions,  generalizing  the  lower  bounds  for  ji.  are  yet  to  be  developed  and  are 
an  important  element  to  be  addressed  in  this  program.  In  the  case  of  the  lower  bounds  for  fi, 
the  available  algorithms  search  for  a  destabilizing  perturbation;  this  would  correspond  in  our 
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case  to  a  value  of  A  such  that  the  previous  kernel  is  nontrivial.  If  such  a  value  is  found  in 
the  unit  ball,  it  would  give  (together  with  the  disturbance  signals  obtained  from  the  kernel) 
the  values  of  uncertainty  which  account  for  the  data. 

This  last  point  has  an  important  implication:  an  extension  of  the  /i  lower  bound  alp- 
rithms  not  only  contributes  to  the  model  validation  problem,  but  also  to  system  identification 
question:  if  some  of  the  uncertainty  in  the  A  structure  is  real  parametric,  it  would  give  the 
values  of  the  parameters  which  are  used  to  fit  the  data.  This  reinforces  the  argument  that 
system  identification  and  model  validation  are  essentially  equivalent  problems. 

5.2.3  Incorporating  time  domain  data  and  dynamical  models 

The  previous  section  was  based  on  static  representations  for  models  and  data,  with  no  explicit 
time  variable.  System  Identification  ,  on  the  other  hand,  deals  with  dynamical  models  and 
observations  of  a  given  variable  across  time.  However,  since  the  horizon  is  finite  it  is  a  concep¬ 
tually  simple  task  to  transform  the  time  axis  into  a  spatial  axis,  by  writing  vectors  with  the 
whole  time  history  of  a  given  variable,  and  large  matrices  which  include  the  equations.  More 
details  of  this  procedure  can  be  found  in  Section  4  and  references  therein.  This,  in  principle 
transforms  the  dynamic  problem  to  the  static  situation  described  before.  However,  computa¬ 
tional  difficulties  arise  from  building  these  large  problems,  and  efficiency  considerations  are 

dominant.  ^  . 

An  alternative  approach  for  LTI  systems  is  to  pose  model  validation  problems  in  the 

frequency  domain;  this  requires  the  availability  of  frequency  domain  data  in  the  form  of 
signal  values  at  a  number  of  frequency  points.  Since  this  is  not  usually  available,  a  previous 
estimation  of  a  number  of  frequency  points  from  the  time-domain  data  is  required,  and  gives 
the  results  only  approximate  value.  From  this  point  onwards  the  problem  can  also  be  turned 
into  a  static  one  in  the  same  way  as  before. 


6  Experim©iital  V^orkt  Active  Control  of  Combustion  Insta¬ 
bilities 

6.1  Introduction 

Oscillations  excited  by  heat  release  mechanisms  inside  combustion  chambers  are  one  of  the 
factors  affecting  the  performance  of  a  diverse  group  of  propulsion  systems.  In  practice,  in¬ 
stabilities  are  avoided  by  changes  in  the  design  of  the  combustion  chamber,  or  by  the  use  of 
passive  stabilization  at  the  expense  of  the  efficiency  of  the  combustion  process. 

Although  the  idea  of  actively  controlling  these  instabilities  has  been  proposed  previously 
([10],  [19],[27]L),  it  has  never  been  implemented  in  a  practical  system.  Recent  research  efforts 
have' applied  the  technique  to  small  scale  experimental  combustors.  However,  this  efforts  are 
all  restricted  to  controlling  one  mode  of  oscillation.  The  successful  implementation  of  active 
instability  suppression  in  cavities  with  complex  geometries  requires  the  ability  to  control 
several  closely  spaced  resonance  modes. 
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We  concentrate  our  effort  on  the  design  of  multi-mode  controllers  for  oscillations  in  a  Rijke 
tube.  To  achieve  this  we  will  use  modern  optimal  controUer  design  techniques  to  modify  the 
linear  acoustic  dynamics  of  the  tube. 

6.2  Experimental  Setup 

The  Rijke  tube  is  one  of  the  simplest  experimental  manifestations  of  thermally  driven  acoustic 
oscillations.  It  consists  of  a  cylindrical  tube  fitted  with  a  heating  element.  The  heater  is  placed 
one  quarter  of  the  way  along  the  tube,  and  its  effective  cross  section  covers  most  of  the  section 
of  the  tube.  The  tube  is  held  horizontally,  and  a  fan  maintains  an  average  air  flow  through 


Fan 


Figure  13:  Description  of  the  Experimental  setup 

A  loudspeaker  is  used  as  actuator.  It  is  placed  with  its  main  plane  parallel  to  the  tube  axis, 
and  its  center  aligned  with  the  intake  opening  of  the  tube.  The  sensor  used  is  a  miniature, 
FET  based  microphone  placed  inside  the  tube.  Both  the  sensor  and  the  actuator  are  placed 
at  the  intake  end  of  the  tube  to  make  controller  design  easier.  Figure  13  gives  a  schematic 
diagram  of  the  experimental  setup. 

Controllers  are  implemented  using  a  real  time  computer.  Analog  to  digital  and  digital  to 
analog  converters  in  the  real  time  system  are  clocked  at  10  kHz,  and  all  sampled  signals  are 
low  pass  filtered  at  3  kHz  before  sampling. 

Both  the  mean  air  speed  and  mean  heat  input  can  be  varied  to  allow  for  controller  robust¬ 
ness  measurements.  Also,  the  tube  has  been  built  with  an  emphasis  on  the  use  of  standard 
parts.  This  makes  results  obtained  easier  to  reproduce  and  verify  and  will  also  help  in  estab¬ 
lishing  a  meaningful  benchmark  for  controller  design  techniques. 

6.3  Proposed  Model 

For  controller  design  purposes  we  assume  that  the  acoustic  dynamics  of  the  tube  are  linear  at 
the  noise  levels  present.  Three  energy  sources  drive  the  pressure  oscillations  inside  the  tube: 
noise  from  the  environment,  the  heat  release  at  the  heating  element  and  the  speaker.  The 
heat  release  is  a  function  of  the  acoustic  speed  at  the  position  of  the  heating  element.  Since 
the  heat  release  does  not  depend  on  the  direction  of  the  air  flow,  this  driving  is  nonlinear. 
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Figure  14:  Proposed  Model 

The  acoustic  dynamics  are  passive  and  therefore  stable,  but  under  certain  conditions 
the  feedback  provided  by  the  heating  element  (thermo-acoustic  feedback)  makes  the  overall 
resulting  system  unstable.  The  inherent  non-linearity  of  this  feedback  is  one  of  the  factors 
that  will  limit  the  amplitude  of  the  oscillations  established. 

The  objective  of  the  controller  is  to  prevent  the  thermo-acoustic  feedback  from  shifting 
the  system  from  a  stable  into  an  unstable  one.  This  objective  is  accomplished  by  modifying 
the  acoustic  behavior  of  the  tube. 

Another  approach  to  this  problem  is  to  act  on  the  heat  releasing  element  and  modify 
the  thermo-acoustic  feedback  characteristics,  however  this  requires  an  accurate  model  of  the 
underlying  phenomena  in  the  heat  exchange  process.  The  mechanical  design  of  the  actuators 
for  this  approach  is  also  harder. 

6.4  System  Identification 

As  has  been  stated  in  the  preceding  section  we  only  need  a  linear  model  for  the  acoustic  dy¬ 
namics  of  the  system.  We  lump  the  dynamics  of  the  amplifier,  speaker,  tube  and  microphone 
into  a  single  system,  and  we  measure  its  transfer  function.  This  is  done  by  applying  a  periodic 
chirp  signal,  with  energy  content  in  the  relevant  band,  to  the  audio  amplifier  input  (d),  and 
measuring  the  output  of  the  microphone  (y).  To  obtain  a  better  signal  to  noise  ratio,  we 
repeat  the  experiment  30  times  and  average  the  results. 


Figure  15:  System  I.D.  and  controller  setup 
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For  computational  purposes  we  approximate  this  transfer  function  with  a  state  space 
realization  of  the  modal  system: 


d{s)  ~  ‘  ;“i  +  ul 

Figure  16  compares  the  measured  transfer  function  with  this  analytical  approximation. 


Figure  16:  Measured  transfer  function  and  analytical  approximation 

In  order  to  carry  out  these  measurements,  we  need  to  place  the  system  at  a  stable  operating 
point.  This  is  best  achieved  by  using  a  higher  air  flow  rate  than  the  one  at  which  the  acoustics 
inside  the  tube  are  unstable.  Figure  17  describes  qualitatively  this  situation.This  implies  that 
the  controller  will  be  used  on  a  system  different  than  the  one  modelled,  and  therefore  places 
an  additional  requirement  on  its  robustness:  the  region  in  which  the  required  performance  is 
achieved  has  to  include  a  neighborhood  of  the  point  o. 


Figure  17:  Qualitative  description  of  robustness  requirements 


6.5  Controller  Design 

As  can  be  seen  in  figure  16  the  transfer  function  presents  resonance  modes  with  a  high 
selectivity  factor  Q  defined  as: 

^  _ _  ^tnax 


22 


where  Umax  is  the  frequency  at  which  the  relative  maximum  is  achieved  and  Aw  is  the  3dB 
bandwidth.  (The  relative  value  of  Q  can  be  estimated  by  the  aspect  of  the  peak:  the  sharper 
the  peak,  the  higher  the  selectivity  factor.) 

In  the  modal  representation  of  the  system,  a  high  selectivity  factor  of  a  resonance  mode 
corresponds  to  a  damping  factor  6  much  smaller  than  one,  or  equivalently  to  poles  close  to 
the  imaginary  axis. 

The  performance  objective  of  the  controller  is  to  simultaneously  lower  the  Q  factor  of  the 
first  two  resonance  modes.  We  follow  the  Loop  Shaping  controller  design  procedure  described 
in  [34]  ,  which  consists  of  3  steps. 


Figure  18:  Loop  Shaping  design  setup 

First,  we  select  the  open  loop  loop-shape  L  needed  to  obtain  our  performance  specification. 
In  our  case  we  require  L  ^  1  over  the  frequency  range  we  want  to  control,  as  this  leads  to 
a  fiat  closed  loop  response  in  that  range.  Additional  factors,  e.g.  the  available  bandwidth 
of  the  controllers,  will  restrict  the  achievable  performance  and  have  to  be  taken  into  account 
while  designing  the  loop-shape.  We  then  design  a  Loop-Shaping  weight  W^s  such  that  L{s) 
Wls{s)H{s). 

The  loop  shaping  weight  we  use  is: 

7 

~  (s+  1500)(5-f  3000)2 

Next,  we  design  a  compensator  to  make  the  closed  loop  system  stable.  The  compensator 
design  must  also  make  sure  that  command  efforts  are  kept  within  operational  limits  of  the  ac¬ 
tuators.  This  is  achieved  by  designing  the  controller  which  stabilizes  the  plant  and  minimizes 
IIGlloo  where  G  is  defined  by  (see  figure  18)  : 


Finally  the  actual  feedback  controller  is  constructed  by  combining  the  Loop  Shaping 
weight  and  the  compensator. 

For  our  application  we  also  reduce  the  order  of  the  controller  using  balanced  model  reduc¬ 
tion,  and  convert  the  continuous  time  controller  into  a  discrete  time  controller.  The  controller 
we  designed  following  this  procedure  has  order  14. 


23 


6.6  Results 

Figures  19  and  20  compare  the  open  and  closed  loop  transfer  functions  of  the  system.  The  Q 
factors  of  the  first  two  modes  were  reduced  from  34  and  42  to  10  and  20  respectively.  This  is 
enough  to  prevent  the  system  from  going  into  limit  cycles.  The  controller  can  also  bring  the 
system  back  from  unstable  to  stable  operation. 


Figure  19:  Comparison  of  open  and  closed  loop  transfer  functions 


Figure  20:  Comparison  of  open  and  closed  loop  transfer  functions  for  the  first  two  modes 

Figure  21  shows  a  power  spectrum  of  the  pressure  waves  inside  the  tube  when  air  flow  and 
heat  input  mean  values  are  set  to  an  unstable  operating  point.  When  the  controller  is  used 
the  limit  cycle  power  levels  are  reduced  by  at  least  60db  in  the  first  two  modes.  In  Figure 
22  we  show  an  interval  of  the  time  record  of  the  output  of  the  microphone  for  both  the  open 
and  closed  loop  configuration. 
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